function W=W(lambda,K_lower,ell,beta,gamma,phi,R,data)

%% comment out if not checking code
% lambda = lambda_vec(3);
% K_lower = 25000;
% beta = beta_vec(3);

%% simulated distribution
Z = [0:500:65000];
% Z = [0:1000:200000];
% polynominal = [1,-.00001,0,0,0,0]; %uniform
polynominal = [1,0,0,0,0,0]; %uniform

weight = zeros(size(Z));
for i=1:length(polynominal)
weight = weight+polynominal(i).*Z.^(i-1);
end
weight(weight<0)=0;
weight = weight./sum(weight);

%% empirical distribution
% data = readtable('covideidlfrequency.csv');
Z = data.bin;
weight = data.cf;

%%
alpha = 1/3;nu=.83;sigma=alpha*nu/(1-(1-alpha)*nu);
func = @(x)((1/sigma*(1-(1-x)^sigma)-x)*R - lambda);

% func = @(x)((1-(1-alpha)*nu)/(alpha*nu)*(1-(1-x)^(alpha*nu/(1-(1-alpha)*nu))-x)*R - lambda);
%func = @(x)((1/alpha*(1-(1-x)^alpha)-x)*R - lambda);
theta = fsolve(func, 0.3);
K_upper = K_lower/(1-theta);

K_star = Z;
K_star(Z>=K_lower&Z<=K_upper) = K_lower;


A = R.*Z.^(1-(alpha*nu/(1-(1-alpha)*nu))).*(alpha*nu/(1-(1-alpha)*nu))^(-1);

W_optimizing = (-phi+beta.*K_star - lambda.*Z).*(Z>K_upper)+A.*K_star.^sigma-ell*K_star;
W_nonopt = (-phi+beta.*Z - lambda.*Z).*(Z>K_lower)+A.*Z.^sigma-ell*Z;

% W_optimizing = (-phi+beta.*K_star - lambda.*Z).*(Z>K_upper)+A.*K_star.^(alpha*nu/(1-(1-alpha)*nu))-ell*K_star;
% W_nonopt = (-phi+beta.*Z - lambda.*Z).*(Z>K_lower)+A.*Z.^(alpha*nu/(1-(1-alpha)*nu))-ell*Z;



%W = sum((W_optimizing*(1-gamma)+W_nonopt*(gamma)).*weight);

 W = sum((W_optimizing*(1-gamma)+W_nonopt*(gamma)).*weight)./sum(Z.*weight);
end

